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^ Abstract 

fr^ \ Space fractional convection diffusion equation describes physical phenomena where 

particles or energy (or other physical quantities) are transferred inside a physical sys- 

^ ' tern due to two processes: convection and superdiffusion. In this paper, we discuss the 

^ \ practical alternating directions implicit method to solve the two-dimensional two-sided 

1-p^ . space fractional convection diffusion equation on a finite domain. We theoretically prove 

and numerically verify that the presented finite difference scheme is unconditionally von 

Neumann stable and second order convergent in both space and time directions. 
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fr^ ■ 1. Introduction 

Relaxing the restriction of the boundedness of the second moments, fractional deriva- 



/\ \ tives naturally appear to characterize anomalous diffusion 13|, usually the time fractional 

H ■ 
_Cy_' derivative is used for describing the subdiffusion and space fractional derivative for the 

superdiffusion. More often, diffusion corresponds to a power law, {x'^(t)) ~ Dt"', where D 
is the diffusion coefficient and t the elapsed time. In a classical diffusion process, a = 1. 
If a < 1, the phenomenon is called subdiffusion. If a > 1, the particles undergo superdif- 
fusion. For a = 1, the position probability density of the particles satisfies the classical 
diffusion equation with first order time derivative and second order space derivative; for 
a < 1, the position probability density of the particles is the solution of the time fractional 
diffusion equation with a-th order time derivative and second order space derivative; for 
a > 1, the position probability density satisfies the space fractional diffusion equation 
with first order time derivative and (3 — a)-th order space derivative. 
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The space fractional advection diffusion equation describes tfie physical phenomena 
involving two physical processes: convection and super diffusion, i.e., for the equation, 
besides the a-th order space fractional derivative term there exists the classical first order 
space derivative. This paper focuses on the numerical algorithm of the following two- 
dimensional two-sided space fractional convection diffusion equation 

OU { X V t] 

^^-^— =d+ (x, y) ^^Dy{x, y, t) + d_{x, y) ^D'^^u{x, y, t) 

+ e+{x,y) y^D'^u{x,y,t) + e_{x,y) yD'^^u{x,y,t) y^-^> 

+ gix, y)u^ix, y, t) + h{x, y)uy{x, y, t), 

where {x,y) eQ = {xl.xr) x {yL.yR), <t <T, the fractional orders 1 < a, (3 < 2; the 
function s{x,y,t) is a source term; and the diffusion coefficients d^{x,y) > 0,(i„(x, y) > 
0,e^{x,y) > and e_(a:, y) > 0. The initial and boundary conditions are, respectively, 
taken as 

u{x,y,0) = uo{x,y) for (x, y) G fi, (1.2) 

and 

u{x,y,t)\Qn = B{x,y,t). (1.3) 

The left and right Riemann-Liouville fractional derivatives of order /i(0<n— l</i<n 



and n is an integer) are, respectively, defined by [4, [ij, ll5| 



m 



XL 



D>(x) = D%,D-^--Mx)l (1.4) 



and 



where 



.D^,^u{x) = E%D-^--^-^u{x)l (1.5) 
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.,D-/u{x) = ^^ I {x- O'-MOd^, ^ > 0, 



XL 



cD->{x) = ^ £' a - xy-MOdi. y > 0, 



and 



d 






d 




E^ 


-D = 




fix' 






dx 



The Griinwald-Letnikov definitions for the left and right fractional derivatives are, respec- 
tively, given as 

^^DM^) = Jim^^5^(-l)^rU(x-^M, (1-6) 



and 

M. 



^dK^) = Jim^ ^ E(-l)^ (fl <x + th.), (1.7) 



where h+ = (x — xl) jM^, h^ = (xr — x)/M_, and M+ and M_ are positive integers. 
The Riemann-Liouville and Griinwald-Letnikov derivatives are equivalent 7| under the 
assumptions that the function performed are sufficiently smooth. 

Using the formulae (11. 6p and (11. 7p to discretize the space fractional derivatives is 
a nature idea for designing the numerical schemes. Unfortunately, this usually leads 
to unconditionally unstable finite difference schemes, but Meerschaert et al successfully 
circumvent this difficulties by modifying the formulae (II. 6p and (II. 7p to obtain the so- 
called shifted Griiwald formulae llO|. Based on these shifted formulae, Meerschaert et al 

id, id, 



did a series of works for numerically solving space fractional diffusion equations 



n 



20| . Sousa proposes another way to approximate the fractional Caputo derivatives |l6l . 
which can obtain second order accuracy. More recently Sousa et al further discuss using 
;he similar idea to discretize Riemann-Liouville fractional derivatives in infinite domain 



18| . By using different ideas, Liu's group discusses the finite difference methods for 
fractional partial differential equations with Riesz space fractional derivatives J2l| and 
the difference methods for the space-time fractional advection-diffusion equation J9|. For 
the discret^ation of the time fractional derivatives, usually the different challenges will 
be met 



|il,3,3l8|. 



Here we will combine the alternating directions implicit (ADI) method 
with Crank-Nicolson scheme to design the finite difference scheme for the two-sided two- 
dimensional space fractional advection diffusion equation (II. ip . The numerical scheme 
will be theoretically proven and numerically verified to be unconditionally von Neumann 
stable and second order convergent. 

The paper is organized as follows. In Section 2, we derive the linear spline approxi- 
mation to the right Riemann-Liouville fractional derivative, and the full discretization of 
(II. ip is presented, where the Crank-Nicolson scheme and the alternating directions im- 
plicit method are combined together. Section 3 does the detailed theoretical analyses for 
the consistency and stability of the given schemes. To show the effectiveness of the algo- 
rithm, we perform the numerical experiments to verify the theoretical results in Section 
4. Finally, we conclude the paper with some remarks in the last section. 

2. Discretization Schemes 

We use three subsections to derive the full discretization of (ll.ip . Since the linear 



spline approximation for the 
terval can be easily got from 



eft Riemann-Liouville fractional derivative in the infinite in- 
18| , we further derive the linear spline approximation for the 
right Riemann-Liouville fractional derivative and make some remarks on the relationship 
between the discretization schemes of left and right Riemann-Liouville fractional deriva- 
tives in the first subsection. Then in the second subsection, we present the scheme for the 



one dimensional case of f ll.ip . The third subsection detailedly provides the full discrete 
scheme of the two-dimensional two-sided space fractional convection diffusion equation. 



2.1. Discretizations for the left and right Riemann-Liouville fractional derivatives 

Let the mesh points Xi = xl + iAx, i = 0,1, . . . , N^, where Ax = {xr — xl)/N,j. is the 
uniform space step. Taking a G (1, 2) in the left Riemann-Liouville fractional derivative 
(11. 4p . its approximation operator S'^^^^u^ has second order accuracy in a bounded domain 
(proved in next section), where m" denotes the approximated value of u{xi,tn) and the 
left fractional approximation operator, which can be obtained by truncating its infinite 



version 



18|, is defined as 



where 



and 



Oi,fc 



8' v"- 



1 



r(4 - a) Ax' 



E 

fc=0 



'^kPi,kJ 



Plk 



0'i-l,k — '^0'i,k + Oj+l,fc5 k < i — 1, 

— iQii + flj_(_i^j, K ^ I, 

0'i+l,i+l) k = I -\- 1, 

0, k>i + l. 



a - 1)3-" - i'^-'^a - 3 + a) 



k = 0, 

{i-k + 1)3-° - 2(i - A;)3-" + {i - k - if-"", l<k<i-l, 
1, k = i. 



(2.1) 



(2.2) 



Here in the finite interval xl < x < xr with a G (1,2), we further derive the linear 
spline approximation for the right Riemann-Liouville fractional derivative defined by 

I q2 fxa 



cD:,u{x,t) 



For a fixed time t, denote 



Ia{x) 



then 



r(2 - a) dx^ 



n(e,t)(e-a;)^-"cie 



r(2 



^ r\(e,t)(e-x)^-"t/e, 

-a) Jx 



M^Mx.t) 



dx'^ 



-'-a \X) 1 



and we can do the following approximation at Xi, 



dx'^ 



-^a\Xi) — 



1 



Ax^ 



(2.3) 



(2.4) 



(2.5) 



[X,(xi„i) - 2X^{xi) + X„(x,+i)], 1 < i < AT, - 1. (2.6) 



For each a;,-, take 



r(2 - a) A^ 



where the sphne Si{^) is defined by 



N^ 



SiiO = ^u{xk,t)s,^kiO^ 



with Sjfc(^), in every subinterval [xk-i,Xk+i], for i + 1 < k < N^ — 1, given as 



■5i,fc(0 = ^ 



•^fc -^fc— 1 






0, 



Xk-l <^<Xk, 



Xk<i< Xk+1, 



otherwise, 



and for k = i and k = N^, Si^kiO taken as 

Xi+l - ^ 



, Xi <^ <x 



SiAO = < 



and 



0, otherwise 



« + l7 



\N^{0 = ' 



Xn^ - Xn^_i ' 



0, 



XTV^-l < ^ < Xtv, 



otherwise. 



According to ( 12. 7p and (12. 8p . we have 

1 



-^a\Xi) — ^a\Xi/ 



r(2-a) 



a^i+i 



M(x„t) / s,MW-^^) di 



N^^l 



V{2-a 



Xk + l 



+ 



k=i+l -'Xk-i 

u{xN^,t) Si^NAi){i-XiY~''di, 



and 






T{2-a 






r(4-a) 



A ™2-a 



xjv^ A 2-a 



r(2 - a) J^^^_^ r(4 - a 



'i,Afr. 



(2.7) 



(2.8) 



+ p.._ ^ 5^ u{xk,t) / s^,fc(0(e-x,)^-"t/e (2.9) 



(2.10) 



where 



bi 



i,k 



Then 



1, k = i, 

(A; - i + 1)3-" - 2(A; - i)3-" + (A; - i - 1)3-", i + 1 < A; < iV^. - 1, (2.11) 

(3 - a - iV, + ^)(iV, - tf~^ + (iV, - t - 1)3-", k = N^. 



-^a\-^i) — -^a\-^i) 



Ax 



2-a 



N^ 



T{A-a 



y^^u{xk,t)bi^k, 



(2.12) 



k=i 



and (12.61) can be written as 






-^a\-^i) 



Ax^ 



[X„(Xi_i) - 2X„(Xi) + Ia(Xi+i)] 

J^ M(Xfc,t)6i_i,fc-2^M(Xfc,t)6j,fc+ ^ M(Xfc,t)6i+i,, 



A^a- 



r(4 - a)Ax 



.k=i-l 



k=i+l 



r(4 - a)Ax^ 



k=i-l 



(2.13) 



where 



k <i-l. 



Qi,k = < 



k = i — 1, 
k = i. 



(2.14) 



0, 

—'2bi^i + bi-i^i, 
^ bi^i^k - '2bi^k + &i+i,fc, i + 1 < A; < A^^. 
Denoting m" as the approximated value of M(xj,t„), we can define the right fractional 
approximation operator as 

E <€k. (2.15) 



A' 7," ■= 



which has second order accuracy for approximating (12.31) (proved in the next section). 
Remark 2.1. Denoting [/" = [m", Wg,--- ^u]^ _-^]^, and rewriting (12.1 p and fl2.15p as 
matrix forms 6'^ , ^f^" = AU"' + 61 and 5^ _2,.f^"' = BU^ + 62, respectively, then there exists 
A = B^. 

2.2. Numerical scheme for one- dimensional fractional convection diffusion equation 

We now examine the full discretization scheme to the one-dimensional two-sided frac- 
tional convection diffusion equation 

du{x,t) 



dt 



d+{x) ^^Dy{x) + d_{x) ^D'^u{x)+g{x)u^{x,t) + s{x,t). (2.16) 



In the time direction, we use the Crank-Nicolson scheme. The central difference formula, 
left fractional approximation operator f l2.ip . and right fractional approximation opera- 
tor (12.151) are respectively used to discretize the classical first order space derivative, left 
Riemann-Liouville fractional derivative, and right Riemann-Liouville fractional derivative. 
Taking the uniform time step At and space step Ax, and setting d+^i = d+{xi),d-^i = 
d_{xi),gi = g{xi), and s""^^^^ = s(xi, t„+i/2), where tn+1/2 = (^n + t„+i)/2, the full dis- 
cretization of f l2.16p has the following form 



At 



r(4-a)Ax^ 



i+l 



^Plkd+,i- 



n+l 



+ ul 






n.+l 



.k=0 



k=i~l 



5i,fc'^-,«" 



nv 



+ 



g, fu:+,' + u-^. 



u 



n+l 

i-1 



(2.17) 



u 



i-l 



2Ax V 2 

Similar to (12. ip and (I2.15p . we define 



^71+1/2 



u'Ui - u 



i-1 






2Aa; 

i+l "'i- l 

2Aa: 
d+,i 



9i\ 



i+l 



^ 7/" •= - 

"'+^ * • r(4-a)Ax^ 



7 ,^kPi,k'^ 



(2.18) 



S" m" ■= ^ 

°-" ' ■ r(4-a)Ax 



fc=0 



a / V k^%. 



fc' 



k=%-\ 



then (I2.17P can be expressed as 

At 



^-T« 



5! 



D" ) 



«"+' 



At 



1 + ^ (^a,+x + (^a.^x + Dl^x) 



(2.19) 
for i = 1,2,..., N^ — 1, associated with the boundary conditions Uq and u^ . Putting 
^i = 2r{4-l)Ax- d+,i, rii = 2r(4^a)Ax" ^-,1, and 7^ = -^gi, the system of equations given by 
(I2.17P takes the form 

(/ - A)t/"+i = (/ + A)f/" + At5"+^/^ (2.20) 



where / is the identity matrix, and 



[/" 



b" n" n" 



,u 



n iT 



^n+1/2 ^ [^' 



n+1/2 n+1/2 n+1/2 



n+l/2iT 



n-f-i/z n-f-i/z-i 

, S2 , • • • , S^ J 



and the matrix entries Aij for 2 = 1,..., N^ — 1 and j = 1, . . . , A'^,. — 1 are defined by 

CiP?,i + ViQ.tv J = h 

iiPli+i + ViQti+i -Ih 3 =i + 1, (2-21) 

6Ki' J < ^ + 1, 



A, 



«j 



and Aqa = 1, Ao,i = for i = 1, . . . , iV^; A^^^n^ = 1 and A^^^i = ior i = 0, . . . ,N^ - 1. 



2.3. ADI scheme for two-dimensional two-sided fractional convection diffusion equation 
Under the direction of discretizing the one- dimensional case of f ll.ip in the last sub- 
section, we use ADI [l9|, y, |5| to numerically solve fll.ip . First we introduce and list the 
denotations (some of them already given above) that will be used in the following: 



D' m". •= 



2Ax ' 






2Ax 



j+i 



1 " " 



N^ 






S" 7/" ■ 



" a,-x"'ij ■ 



"+,i I 



E 

k=i~l 
i+1 



'^k,jQi,k'i 






r(4 - a)Ax 



fc=0 

A'. 



7 , '^k,j%ki 



n 7/" •- 



'^ij+l ^i,i-l . 



2Ay 



7-)/' „,n ._ ^».J+1 "^.i-U 
^P,y%j ■— oA„. "m' 



2A2/ 



J+1 



1 

^P,+y^i,j '■= r(4 _ «) Aw" ^ ^i,kPj,k'^ 
y 7/^^ — ^ V 7/" «^ ■ 



r(4-/3)Ay° 



i+1 



A'' ,," ._ ''+.M \^.," „^ 



r(4-/3)Ay^ 



fcl 



fc=0 



-,t,J 
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k=j-l 



(2.22) 



Analogously we still use the Crank-Nicolson scheme to do the discretization in time di- 
rection. Taking u"^- as the approximated value of u{xi,yj,tn), d^^ij = d+{xi,yj), d^^ij = 

d-yXiyHj), ^+,i,j — ^+[^i^yj)y ^-,i,j — ^-[^i^Uj): 9i,j = dy^iyVj): i^ij = h[Xi,yjj, tn+l/2 = 

(t„ -F t„+i)/2, s'^j = s{xi, yj, tn+1/2), Ax = {xr - xl)/N^, and Ay = {yn - yL)/Ny, for 
the uniform space steps Ax, Ay and time step At, the resulting discretization of (II. ip can 
be written as 






1 



At 



r(4-a)Ax^ 
1 






,«j 



.fc=0 



+ 



+ 



T{A-P)AyP 



/ . Pj,k'^+,i,j 
.k=0 



k=i-l 






2 

Kt' + Kk 



r, I 7,"+l -L 7," 7,"+l -I- 7," 



2Ax 



+ 



^*J I ^«J+1 ~'~ "^i.i+l ""ij-l ~'~ ""ij-l 



+ S 



2A2/ 

n+1/2 



(2.23) 



Using the notations (12.221) . we have 



At 



{K,+x + K,-x + D'^,^) 



At 



5'' 



l3,+y ~T- ^I3^_y 



6' 



D'' 



I3,y 



u. 



n.+ l 
i,3 



1 + ^ {^'L,+x + ^'L,_x + D'^,x) + -TT {^'fi,+y + ^'fi,_y + D'^^y) 



U 



«J ' *J 



."f/'At. 



Further define 



(^a.x := ^^ , .^ + ^^ .^ + -0" 



a,4-X ' "a,^x 



a,xJ 



Sl3,y ■= S" , „ + 5'^ „ + D'^ 



(2.24) 



(2.25) 



thus, (12.241) may be rewritten as 



1 7rOa,x —0^ 



At 



2 ~(^'yj ^id 
At 



n+l 



1 + ^6^,x + -Yh,y ] <, + s^r ^- 



^n+l/2 



(2.26) 



For the two-dimensional two-sided fractional convection diffusion equation (II. ip . the rel- 
evant perturbation of (I2.26P is of the form 



1-^A 



1 + — 5, 



At 



1 - ^5p,y ] <:t 



1 + ^h. ) <. + <r''A^- 



(2.27) 



The system of equations defined by (I2.27P may be solved by the following ADI (Peaceman- 
Rachford type) scheme 191 ]: 



At 



1 



2 

At 



At 

y 



5a,x ) Ulj = ( 1 + —5/3,2; ) K,j + ^^^'i ' 



At 



,"+1/2. 



5/3,2^ ^ 



,n+l 



1 ^^r 

1 + — 5n 



2 ^^'y / "«J \ " ' 2 '""''' y ^*'-J' ~^ 2 '^*'-' 



At 



n+l/2 



(2.28) 



(2.29) 



where m* • is an intermediate solution. For maintaining the consistency, we need to care- 
fully specify the boundary conditions of uly Subtracting (I2.29P from (I2.28p . we obtain 



2<, 



At 



h,v <r + 1 + ^^^. 



At 



y I ^ij' 



(2.30) 



then the boundary conditions for u*j ( i 
given as 



and i = Nx with j = 1, . . . , A'j^ — 1) can be 



u. 



o,j 



u 



N,,j 



1 

2 
1 
2 
1 
2 
1 
2 






n+l 



-(5fl „ 1 Un 



Y^f^^y ) ^O'i 



»i+i 



At 
At 



J/3,2/ 






At 






-f^ 



z?.?/ ] ^S 






(2.31) 



The corresponding procedure is executed as follows: 
(1) First for every fixed y = yk {k = 1, . . . ,Ny — 1), solving a set of A''^. — 1 equations 



defined by (12.281) at the mesh points Xi,i = 1, 



A^:, - 1, to get u* 



fci 



(2) Next alternating the spatial direction, and for each fixed x = Xk {k = 1, . . 
solving a set of A'j^ — 1 equations defined by (I2.29P at the points yj,j = 1, . 
to obtain m^j^- 



.,Ny-l, 



3. Convergence and Stability Analysis 

We show the convergence for one-dimensional and two-dimensional two-sided fractional 
convection diffusion equation by proving the consistency and stability (according to Lax's 
equivalence theorem). 

3.1. Convergence and stability for one dimensional two-sided fractional convection diffu- 
sion equation 
Lemma 3.1. Let u G (7^(0) and ^ G [xfc,Xfc+i], then we have 

,(Ax)2 e{i - e){i - 2e){Axf 



u{0 - s.iO = oa 



-u"{0 + 



3! 



-u 



'(0 + O ((Ax)^) 



where Ax = Xk+i — Xk, = {^ — Xk)/Ax. 



Proof. For any ^ G [xk,Xk+i], using Taylor series expansion at ^, it is easy to get the 
results. 

D 

Lemma 3.2. j^, |lj] Let D^u{x) be continuous in the interval [xL,xji] and < n — 1 < 
fi < n, then for xl < x < Xr the following holds 



.,D>(x) = D%^D-^'^-^Mx)] =.p-^'--'\DMx)] + J2 



'x — XlY ^ 



fc=0 



r(-/i + A:+l) 



D'f{xL). 
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Lemma 3.3. Let D^u{x) be continuous in the interval [xi,Xij] and 0<?2 — l</i<n, 
then for xl < x < xr the following holds 



n~l , \k—n 

Xr — X) ^ 



^D^.A^) = E%D-tMx)] =.D-t'^[EMx)] + E f7r!—T-TT^V(x/j). 



fc=0 



r(-/x + fc + i)' 



Proof. Since xD^^u{x) = p7^ J^^ {^ — x)" ^u{^)dC,, i^ > 0, it can be obtained 

1 r^ dif - rY~'' 1 /"^^ 



u-l 



^DZ[Du{x)] - ^^^^/(^ij) 



then 

Viu 



ElD7Zu{^)] =.DZ[Eu{x)] + ^^^V7;^/(^«)- (^.l) 



Using E to differentiate both sides of (13. ip leads to 

E%DZu{x)] = E{ElDZu{x)]} = E{^DZ[Eu{x)]} + ^^Y^^^H^r)- (3-2) 
Replacing u with Eu in fl3.ip . then (13. 2p can be rewritten as 

Repeated iterations establish the desired result. 

D 

Theorem 3.4. Let u E C\n) and xD^^^u G ^^(0), then ^D'^^u{xi,tn) = 5'^,_xUi + 
C((Ax)2), where a G (1,2). 

Proof. In the following we use u{xi) to denote u{xi,tn)- Since i„ is fixed and xDf.~^u G 
C"^(r2), then the following holds 

92 1 

^2:,(x,) = -^^ 

Let e(xj) be the error satisfying 

Xa{Xi-i) - 2Xa{Xi) +Xa{Xi+i) = Ia{Xi-i) - 2Ia{Xi) + /^(Xi+i) + e{Xi), 

then there exists 

xD:^u{x,) = ^^ [/,(x,_i) - 24(x,) + /„(x,+i)] + ^^e(x,) + O ((Ax)2) , 

11 



xD^uixi) = —^Ia{xi) = [I^{xi^i) - 2X^{xi) +X<^(xi+i)] + O ((Ax)^) . 



I.e., 



D"u(x^) = S' ,X + 



1 



x-^x 



a, — xt 



Ax,) + O ((Ax)2) . 



(Ax)2' 

Next we calculate the error e{xi). Extending the definition of u{x) from [x^jXr] to 
[xl,xr + Ax] with the extended part of u{x) being the one given in Theorem A of Ap- 
pendix, then 



Jxi-l 



xr 



Jxi JXi + l 

\u{^ - Ax) - 5,_i(e - Ax))(e - x,)'^"d^ 



xr 



Hi)-s^{mi-x^y-''di 



+ 

I X 



xr-Ax 



(n(e + Ax) - ^,+i(e + Ax))(e - x,)^-"rfe 



2-'fe+i 



J2 I («(e - Ax) - ^,_i(e - Ax))(e - x,)'--d^ 

k=i -^^k 

- 2 Y, r^\u{0 - S,{0){^ - x,)'-^d^ 
k=i -^^fe 

+ £ r''''(M(e + Ax)-5,+i(e + Ax))(e-x,)i-°rfe + c(x, 
k=i "^^fc 



where 
C(x,) 



(«(0-5._i(0)(e-a:.-i)^-"rfe 



XH+Aa; 



(«(0-5.+i(0)(e-a;.+i)'-"rfe 



and C{xi) = 0. Using Taylor series expansion and denoting 6 = (^ — Xfc)/Ax, from Lemma 
3.1 we get 



N^-l 



eiXi 



l-a 



k=i >^^fc 



^(l-^)(l-2^)(Ax) 
3! 



^(l-^)(Ax) 
2! 



K(e-Ax)-2u"(0 + n"(e + Ax)) 



M 



'(e - Ax) - 2u"'{i) + M'"(e + Ax)) + O ((Ax)^) 



di. 



Further using the first mean value theorem for integration and Taylor series expansion, 



12 



~k 



there exist fj , r]^, 772, t]^, and r]^ such that 



eiXi 






^(l-e)(Aa;) 
2! 



■(M"(r/'= - Ax) - 2u"{fi'') + M"(f/'= + Ax)) + 



9{l-9){l-29) 
3! 



x(Ax)3 {u"\fi'' - Ax) - 2u"'{t) + v:'\t + Aa;)) + C ((Ax)"^)] r\ - Xi)^-"^^ 



E 

k=i 



el^l - ff,(A.l ^ (M^,„„„(,,j) ^ „„„,^.)) ^ ^(l^-'Xl - 2C)(Ax)* 



2! 



xK'"(r/3^) - «""(r^4^)) + O ((Ax)^)] H" (^ - x.^-'^d^ 

Jxk ) 

l-^)(Ax)^ ^(l-^)(l-2^)(Ax)^ 



3! 



A^:.-! 



<i"""(')o)i r' ^^"~''' + "■• '""3 """-"' + o ((Ai)') 



a;fc+i 



k=i -^^k 
1 , 



2-a 



^X/j Xj 



\2-a 



where \u""{riQ)\ = max |n""(?7^')|. Therefore e(xi)/(Ax)^ = C((Ax)^), the de- 

i<k<Nx-l,l<j<4: ■' 

sired result is proved. 

D 

Remark 3.1. From Lemma 3.3, it can be noted that u G C^ipi) and xD^~'^u G C^{VL) if 
and only if n G (7^(0) and u^^\xji) = 0, A; = 0, 1, 2, 3. At the same time, if we do the zero 
extension to u at [xR,Xij + Ax], then u G C'^[xL,Xij + Ax] not belongs to C"^[xL,Xij + Ax]. 

Using similar idea we can prove 
Lemma 3.5. Let u G C\n) and x^D^^'^u G C^i^), then x^Dy{xi,tn) = 5l,^^< + 
0((Ax)2), where a G (1,2). 

Remark 3.2. From Lemma 3.4, it can be noted that u G C'^{0.) and x^D^^'^u G C^{0.) 
if and only if n G C^(fi) and u('=)(xl) = 0, A; = 0, 1, 2, 3. 

Example 1. To numerically verify the truncation error given in Theorem 3.4 in a bounded 
domain, consider the function u{x) = sm((l — x)^), x G f2 = (0, 1); and it is easy to check 
that u^^\x)\dn 7^ 0. Denoting /(x) = 2,.D"^m(x), by the algorithm given in 1|, we can 
numerically obtain the value of /(x) at anywhere of the considered rectangle domain with 
any desired accuracy. 

Example 2. To numerically verify the truncation error given in Lemma 3.5 in a bounded 
domain, consider the function u{x) = sm(x^), x G f2 = (0,1); and it is easy to check 
that u^^''{x)\dn 7^ 0. Denoting F{x) = xlD'^u{x), again by the algorithm given in [l|, we 
can numerically obtain the value of F{x) at anywhere of the considered rectangle domain 
with any desired accuracy. 
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Table 1: The inaxiinuin errors (14.11) and convergent orders for the scheme fl2.15p . 



Ax 


a = 1.1 


Rate 


a = 1.5 


Rate 


a = 1.9 


Rate 


1/50 


1.0303e-002 




2.4122e-002 




3.8358e-002 




1/100 


2.7832e-003 


1.8882 


6.4914e-003 


1.8937 


1.03226-002 


1.8938 


1/200 


7.2217e-004 


1.9464 


1.6846e-003 


1.9461 


2.7038e-003 


1.9327 


1/400 


1.8379e-004 


1.9743 


4.2937e-004 


1.9721 


6.9889e-004 


1.9518 



Table 2: The maximum errors (14. ip and convergent orders for the scheme ( 12. ip . 



Ax 


a = 1.1 


Rate 


a = 1.5 


Rate 


a = 1.9 


Rate 


1/50 


1.0303e-002 




2.41226-002 




3.83586-002 




1/100 


2.78326-003 


1.8882 


6.49146-003 


1.8937 


1.03226-002 


1.8938 


1/200 


7.22176-004 


1.9464 


1.68466-003 


1.9461 


2.70386-003 


1.9327 


1/400 


1.83796-004 


1.9743 


4.29376-004 


1.9721 


6.98896-004 


1.9518 



Table 1 and Table 2 numerically verify Theorem 3.4 and Lemma 3.5, respectively, 
and they show that the truncation errors are second order. In fact, watch carefully the 
numerical results, which also confirm Remark 2.1 in some sense. 

Lemma 3.6. The coefficients gf;, defined in (12.141) satisfy Yl ^Zi^fc < 0, 1 < i < A^z — 1- 

k=i-l 

Proof. Taking u{x,t) = l, from ([221 l232l [27[3D . we have 

k=i—l k=i—l 

since u{x,t) = 1 implies its linear interpolation function Si{x) = u{x,t) = 1. 
Then 

/„(Xi_i) - 2Ia{Xi) + Ia{Xi+i) 



r(2-a) 
1 



J Xi^l J Xi J Xi-f-l 

{{XR - x._i)'"" - 2{xR - Xif-'^ + {XR - x,+i)2-") 



r(3-a 

— ^ ((a + Ax)^"" - 2a^-" + (a - ^xf^^) 
r(3 — a) ' 

^^^^ [((a + Ax)^-" - a'-) - (a^-- - (a - Ax)^"")] 

x^-^dx- / x^-'^dx 



(2 - a)r(3 - a) \J^ J^_^^ 
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where a = xr — Xi > 0; obviously, x^ " is a decreasing function, so /q,(xj_i) — 2/q,(xj) + 

laiXi+i) < 0. D 



Using similar method, we can prove 



i+l 



Lemma 3.7. The coefficients p^^ defined in fl2.2p satisfy ^ p"^ < 0, I < i < N^ — 1. 



fc=0 



Theorem 3.8. The scheme (12.191) of the fractional convection diffusion equation fl2.16p 
with constant coefficients and 1 < a < 2 is unconditionally von Neumann stable. 

Proof. Let m" be the approximate solution of m", which is the exact solution of the scheme 
(^J^. Setting e] = uj - u], 1 < j < N^ - I, then from dZTg]) we get the following 
perturbation equation 



1 - Y<^ ^T' = (l + Y^a,. ] el 



(3.4) 



with zero boundary conditions, i.e., Eq = e^^ = 0. Then we can use the Von Neumann 



analysis or Fourier method [3|] to do the stability analysis. Putting Q = ■^tu^^Tax"'^ 



+J' 



Vj = 2f{4^)A^^-J' ^- = ^9j, and assuming 



then (13.41) leads to 



e^+Uu) 



i+i 



/2tt 



N^ 



=£"(u;)e*""^ 



1 - Z1^M^^'"^^^~'^^" + ^^' J2 Qtke""^'"'^^" + 7j (e'"^" - e"'"^") 



k=0 



k=j-l 



e^(uj) 



j+i 



A^. 



1 + ^i Xl^i.fc^''^^'"^'^'^'' + ^^ Yl qlke'''^^-^'^^'' + 7j (e*""^^' - e"^'^^^') 



fc=0 



The amplification factor is 

i+i 



k=j-i 



N^ 



fc=o fc=i-i 

fe=0 A;=i-1 



Next we prove that |Q(w)| < 1, it means that the real part of the element in the big 
parentheses of above equation is negative, that is 

y^p'^i^cos{u{k — j)Ax) < and >^ q",^cos{uj{k — j)Ax) < 0. 

fc=0 k=j-l 
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We can write 

i+i i-2 

"^pjkcosiuik - j)Ax) = {pIj_^ + p°j.+i)cos(wAx) + pjj + ^p°fcCos(w(/i; - j)Ax), 

k=0 k=0 

J2 q-,kCOsHk - j)Ax) = {ql^_, + ql^^,)cos{uAx) + ql^ + J^ g,^,co.(a;(/c - j)Ax). 

k=j-l k=j+2 

From Theorem B of Appendix, we know p'jj^i + P'jj+i > 0, p";, > for /c < j — 2, and 
gj^,_i + qf,,+i > 0, g^^, > for A; > J + 2, then 

i+i i-2 i+i 

Y,PtkCos{u;{k-j)Ax) < (p-^_, +p-^.^,) +p-^. + ^p-, < J2p-^ < 0, 

A;=0 k=0 k=0 

J2 qlkCOs{uj{k - j) Ax) < (g°^.„i + ql^^,) + g°^. + J^ g^^ < $^ g",. < 0, 

k=j-l k=j+2 k=j~l 

where Lemma 3.6 and 3.7 are used. D 

3.2. Convergence and stability for two-dimensional two-sided fractional convection diffu- 
sion equation 

In this subsection, we prove the consistency and stabihty of the scheme (I2.27p -( l2.3ip 
which combine the ahernating directions imphcit scheme with Crank-Nicolson scheme 
together (ADI-CN). According to the Lax's equivalence theorem, the convergence of the 
ADI-CN scheme is naturally obtained. 

Theorem 3.9. The truncation error of the scheme ( l2l27|) is 0{{Ax)^) + 0{{Ayf) + 
0{{AtY). 

Proof. Let u{x,y,t) be the exact solution of the two-dimensional two-sided fractional 
convection diffusion equation fILip . then for the scheme f l2.26p we have 

u^^^ — -u" ■ 

Mn-D' +h-D' \ !f!J_LZM _ e"+i/2 

_du{Xi,yj,tn+l/2, 



o, d+,i,j XLD^U[Xi,yj,tn-i-l/2) — d-^ij xD^^u[Xi,yj,tn+l/2) 

~ ^+,i,j yLDyU{Xi,yj,tn+l/2) — G-,i,j yDy^u{Xi, yj,tn+l/2) ~ 9i,j'^x{Xi,yj,tn+l/2) 

- hijUyixi, yj, tn+1/2) - s{xi, yj, tn+1/2) + ^(Aa;)^ + OiAyf + C(At)l 
The scheme f l2.27p differs from f l2.26p by a perturbation equals to [l9| 
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which may be deduced by distributing the operator products in ( I2.27p . Since (u^j^—u^j) is 
an 0{At) term, it follows that the perturbation contributes an (9((At)^) error component 
to the truncation error of (12.261) . Thus, the scheme f l2.27p has a truncation error also 
0{{Axf) + 0{{Ayf) + 0((At)2). D 

Theorem 3.10. The ADI-CN scheme, defined by (12.261) with constant coefficients, is 
unconditionally von Neumann stable for 1 < a, /3 < 2. 

Proof. Let u"^ be the approximate solution of u"'^, which is the exact solution of the 
scheme ([22SD- Setting el^ = v^ - u],^, I < j < N^ - I, and 1 < m < Ny - 1, then 
from (I2.26P we get the following perturbation equation 

At At 

— a,x- — 

Similarity, putting ^i = 2r(4^)A^^+.i' Vj = 2r(4^)A^^-J' ^J = ^9j and ^^ 

2r(4— /3)Aj/'3 6+,nu Vm ^ 2r(4— /3)Aj/^ ^-,m; Im ^ 4Ay "" ^SSUmmg 

1 



1 —^a,x —^I3,y ] ^i,m 



^ At^ At^ , „ 



(3.5) 



then ( 13. 5p leads to 



£"+i(a;) 



i+i 



1-K.E^M 



'j.m 



^jaj(A;— j)Ax 



/27r 



Afx 



=£ we 



ioiAz ^,~iu>Ax\ 



m E C/^e-^'^-^)^- + 7, (e— ^ - e 



fc=0 
m+l 



fc=i-l 



Af« 



e^E^l.-^^^'=-"^^^^+^~"^ E 



/3 iu}{k—m)A.y 



e-iu) 



fc=0 k=m—l 

A;=0 fc=J-l 



,iivAy _ -iwAy\ 



+ 7j (e*'^^^ - e-^"^^^) 



m+l 



AT^ 



k=m—l 



\ fc=0 

The amplification factor is 



V fc=o fc=i-i 

m+l ^v 



iujAx — iwAz 



) 



^m _^ Pm,fc^ 



/3 ia;(fc— m)Ai; 



+ ^m E ^fce-^'^"™^^^ + 7m (e 



iwAi; —iuAy 



) 



fc=0 



fc=m— 1 

1 + ( 0- E^i.fc^''^^^'^^'^'' + ^i E (llk^""'^^'^^^'" + 7j (e'"'^'' - e"^'^^^) 

fc=i-i 



fc=0 



m+l 



Af„ 



.f \^ r)^ ia;(fc-m)Aj/ 
Sm / ^ PmM^ 



+ V~m E ^™,fc^ 



j8 iu}{k—m)Ay 



+ 7m ( 



giwAy _ g-^wAy^ 



fc=0 



fc=m,— 1 
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Next we prove that |Q(w)| < 1, it means that the real part of the element in the big 
parentheses of above equation is negative, that is 

2_]pJ^k'^os{uj{k — j)Ax) < and >^ q^^.cos{(jj{k — j)Ax) < 0, 

fc=0 k=j-l 

m+l Ny 

Yl Pm,fc^^^(^(^ - m) Ay) < and ^ q'^i^cos{uj{k - m) Ay) < 0. 

k=0 k=m—l 



Similar to Theorem 3.8, we get that the scheme (I2.26p is unconditionally von Neumann 



stable. n 

4. Numerical Results 

Here we verify the above theoretical results including convergent order and stability. 
Introducing the vectors U = [uh{xo, t), . . . , Uh{xiy, t)]'^, where U is the approximated value, 
and u = [u{xQ,t), . . . ,u{xi\i,t)f^, where u is the exact value, in the following numerical 
examples the errors are measured by 

\\U{Ax)-u{Ax)\U (4.1) 

where || ■ ||oo is the /qo norm. 

4.1. Numerical results for one- dimensional two-sided fractional convection diffusion equa- 
tion 

Let us consider the one-dimensional fractional convection diffusion equation fl2.16p . 
where < x < 1 and < t < 1, with the coefficient functions 

(i+(x) = r(3 — a)x", (i_(x) = r(3 — a;)(2 — x)", and g{x) = -x. 

Take the exact solution of the equation as u{x,t) = e~*sm((2x)^)sm((2 — 2x)^), then the 
corresponding initial and boundary conditions are, respectively, u{x, 0) = sin{{2x)^)sin{{2- 
2x)^) and u{0,t) = u{l,t) = 0; and the forcing function 

s{x, t)=- e-*sm((2x)^)sm((2 - 2x)^) - e'^d+^x) ^^D^sin{{2x)^)sin{{2 - 2x)^) 

- d_{x) ,^D2^sm{{2xf)sin{{2 - 2x)^) 

- 64:e-^g{x){x^cos{{2x)^)sin{{2 - 2a;)^) - (1 - xfsin{{2x)'^)cos{{2 - 2xf)), 

by the algorithm given in [l|, we can numerically obtain the value of s(x, t) at anywhere 
of the considered rectangle domain with any desired accuracy. 

In Table 3, we show the scheme (12.171) is second order convergent in both space and 
time directions. 
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Table 3: The inaxiinuin errors (14. ip and convergent orders for the scheme fl2.17p of the one- 
dimensional two-side fractional convection diffusion equation f l2.16p at t=l and At = Ax. 



At, Ax 


Q= 1.1 


Rate 


a = 1.5 


Rate 


a = 1.9 


Rate 


1/50 


2.1180e-003 




1.9815e-003 




1.3809e-003 




1/100 


5.2688e-004 


2.0072 


5.0092e-004 


1.9839 


3.5593e-004 


1.9559 


1/200 


1.3174e-004 


1.9997 


1.2649e-004 


1.9856 


9.1681e-005 


1.9569 


1/400 


3.2913e-005 


2.0010 


3.1851e-005 


1.9896 


2.3523e-005 


1.9625 



4-2. Numerical results for two-dimensional two-sided fractional convection diffusion equa- 
tion 
Consider the two-dimensional two-sided fractional convection diffusion equation (II. ip . 
on a finite domain 0<x<l,0<i/<l,0<t<l, and with the coefficients 

d+{x,y) = r{3-a)x'', t/_(x,|/) = r(3 - a)(2 - x)", g{x,y) = -x, 

e^{x,y) = r(3 - /?)/, e4x,y) = r(3 - /3)(2 - yf, h{x,y) = ^y, 

and the initial condition u{x, y, 0) = sm((2x)^)sm((2 — 2x)'^)sin{{2y)'^)sin{{2 — 2y)^) and 
the Dirichlet boundary conditions on the rectangle in the form u{0,y,t) = u{x,0,t) = 
and u{l,y,t) = m(x, l,t) = for all t > 0. The exact solution to this two-dimensional 
two-sided fractional convection diffusion equation is 

u{x,y,t) = e-*sm((2x)^)sm((2 - 2xY)sin{{2yf)sin{{2 - 2yY). 

By the algorithm given in [l| and above conditions, it is easy to obtain the forcing function 
s{x,y,t) at anywhere of the considered rectangle domain with any desired accuracy. 



Table 4: The maximum errors (14. ip and convergent orders for the scheme fl2.27p - fl2.3ip of 
the two-dimensional two-sided fractional convection diffusion equation f ll.ip at t = 1 and 

At = Ax = Ay. 



At, Ax, Ay 


a= 1.1,13 = 1.1 


Rate 


Q = 1.6,13 = 1.4 


Rate 


Q = 1.9,13 = 1.9 


Rate 


1/25 


9.5946e-003 




8.5313e-003 




1.0232e-002 




1/50 


2.3956e-003 


2.0018 


2.1729e-003 


1.9731 


2.6207e-003 


1.9650 


1/100 


5.9582e-004 


2.0075 


5.5244e-004 


1.9757 


6.6155e-004 


1.9860 


1/200 


1.4915e-004 


1.9981 


1.3959e-004 


1.9847 


1.6796e-004 


1.9778 



Table 4 also shows the maximum error, at time t = 1 and At = Ax = Ay, between the 
exact analytical value and the numerical value obtained by applying the ADI-CN scheme 
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fl2.27p -( l2.3ip . and the scheme is second order convergent and this is in agreement with 
the order of the truncation error. 

4.3. Numerical results for two-dimensional one-sided fractional convection diffusion equa- 
tion 
Considering the two-dimensional two-sided fractional convection diffusion equation 
(11.11) . and taking the coefficients functions as 

d+{x,y) = l, d_{x,y) = 0, g{x,y) = l, 
e+ix, y) = 1, e_(x, y) = 0, h{x, y) = 1, 
then it becomes the two-dimensional one-sided fractional convection diffusion equation, 

du{x,y,t) ^ d"u{x,y,t) d'^u{x,y,t) du{x,y,t) du{x,y,t) ,^^. 

dt dx"' dy^ dx dy 

where 0<x<l, 0<y<l, 0<t<l, and the initial condition u{x,y,0) = 
sin{x^) sin{y^) and the Dirichlet boundary conditions on the rectangle in the simple form 
u{0,y,t) = u{x,0,t) = 0,u{l,y,t) = e~^sin{l)sin{y'^) and u{x,l,t) = e~*sm(l)sm(x^) 
for all t > 0. The exact value to this two-dimensional one-sided fractional convection 
diffusion equation is 

u{x,y,t) = e~^sin{x^)sin{y^). 

By the algorithm given in p| and above conditions, it is easy to obtain the forcing function 
s{x,y,t). at anywhere of the considered rectangle domain with any desired accuracy. 

Table 5: The maximum errors ( 14.1 p and convergent orders for the scheme (I2.27p - (l2.3ip of 
the two-dimensional one-sided fractional convection diffusion equation (14. 2 p at t = 1 and 

At = Ax = Ay. 



At, Ax, Ay 


a = 1.1,^ = 1.1 


Rate 


a = 1.6,/? = 1.4 


Rate 


a = 1.9,13 = 1.9 


Rate 


1/25 


1.1435e-003 




5.0896e-004 




2.63816-004 




1/50 


2.8953e-004 


1.9817 


1.3592e-004 


1.9048 


6.93906-005 


1.9267 


1/100 


6.8091e-005 


2.0882 


3.4877e-005 


1.9624 


1.80646-005 


1.9416 


1/200 


1.5950e-005 


2.0939 


8.8502e-006 


1.9785 


4.67286-006 


1.9507 



Table 5 shows the maximum error, at time t = 1 and At = Ax = Ay, between the 
exact analytical value and the numerical value obtained by applying the ADI-CN scheme 
(I2.27p - (l2.3ip . and the scheme is second order convergent and it corresponds to the order 
of the truncation error. 
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5. Conclusions 

This work provides tlie second-order efficient numerical scheme for the two-dimensional 
two-sided fractional advection diffusion equation on a finite domain. Both the convergent 
order and numerical stability of the scheme are theoretically proved and numerically 
verified. This paper can be considered as the sequel of the works [lO, UjJ, ll2|, uM, 113, llSl, 
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Appendix 

Theorem A. For any given u(x) G C^[xl,xfi], there exists an extension of u{x) defined 
on [xr, xr + Ax], such that the extended u{x) G C^[xli xr + Ax] and 

/ («(e)-s._i(o)(e-x.-i)^-"rfe- / («(o-5.+i(0)(e-x.+i)'-"^e = o, (5.1) 

Jxn-Ax Jxa 

where S'j_i(x) is the linear interpolation function oiu{x) on the interval [xr — Ax, xr\, and 
S'j+i(x) the linear interpolation function of the extended u{x) on the interval [xr, xr+Ax]. 

Proof. Denoting a = XR,h = xr + Ax, and di = £l^_/^^{u{0 - Si-i{0)i^ - Xi-iY~°'d^, 
then the equality (15. ip means to find the extension of u{x) on the interval [a, b] such 
that J^{u{^) — SiJ^i{C,)){^ — Xi+iY'^'dC, = di. Taking u{b) = u{a), then 5^+1(0;) = u{a) 
on the interval [a,b] since Si^i{x) is the linear interpolation of u{x). Now we need to 
prove that there exists V{x) = u{x) — SiJ^i{x) on the interval [a, 6] such that V{a) = 0, 
V{b) = 0, V'{a+) = u'{a-) = c/2, V"{a+) = u"{a-) = d^, V"'{a+) = u"'{a-) = d^, and 
V""{a+) = u""{a-) = 4- 

Suppose that V{x) is composed by {l,x,x'^,x'^,x^,x^,x^,x^}, i.e., 

V{x) = Co + cix + C2X^ -l- csx^ -|- c/^x^ -|- csx^ -|- cex^. 

Using the first mean value theorem for integration, there exists i] G (a, 6), such that 

J a 
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and it can be rewritten as 



Co + ci?7 + C2?7^ + C3?7'^ + C4?7^ + c^rf' + c^y^f 



b X -1 

l-a. 



(e-x,+i)^-X rfi=:rfi 



(again denoted by di, and it is easy to check that /^(^ — Xi+iY~°'d^ 7^ 0). Combining 
all the above requirements, we obtain the linear algebraic equations Ac = d, where c = 
[co, ci, C2, C3, C4, C5, ce]'^, d = [0, 0, di, d2, d^, d4^, d^]^ , and the determinant of the coefficient 
matrix A is 



det{A) 



b 

V 



1 
1 
1 





a 
a 



a 
b 

V 
1 









62 

n2 



a" 

^2 

2a 
2 





3a2 
6a 

6 





a" 

4 
rj 

Aa^ 

12a2 

24a 

24 



a^ 
b' 
rf 

20a3 
60a2 
120a 



?7" 

Ga^ 

SOa"^ 

120a3 

360a2 



1 






Using the formula a" — 6" 



rj — a 

2a 

2 





(a — 6) (a 



n-l 



63 -a3 
r/^ — a^ 

6a 
6 


6a"-2 



4 4 
77 — a 

4a3 

12a2 

24a 

24 

52^n-3 _ 



a 

65-a5 

?7^ - a^ 

5a^ 

20a3 

60a2 

120a 

. + ¥'-' 



a 



a 

b'^-a^ 

rf -a^ 

6a^ 

30a^ 

120a3 

360a2 

6""^), taking 



out the common factors r] — a and b — a, and repeating the steps, finally we get 



det{A) =288{r]-af{b 



I a a a^ 

1 2a 3a2 

1 3a 

1 







--288{r]-b){7]-a)\b-af. 



a" 
4a3 
Qa' 
Aa 

1 







a" 

Sa^ 

lOa^ 

lOa^ 

5a 

1 





a" 
6a^ 
15a^ 
20a3 
15a2 
6 + 5a 
rj — b 



Since rj ^ a, rj ^ b, and a 7^ 6, we get det{A) ^ 0. Then with the given basis functions {1, 
X, a;2, x', x^, x'^, x^, x^}, there exists a unique extended part of u{x), which is V{x) +u{a) 
with X G [xr, xr + Ax]. The extended u{x) G C^fx^, xr + Ax] and (15. ip holds. D 
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Theorem B. For the coefficients gf^. defined in fl2.14p and p";. defined by (12.21) . the 
following hold: g°^ > ioi k > i + 2, p"^ > ior k < i — 2, and p'^^^i = q^i_i = 1, 



Pti = €i 



-4 + 23-", p-^_^ = g° _^, = 6 - 25- + 33-°. 



Proof. Similar to the proof given in 18| for the left Riemann-Liouville fractional deriva- 
tive, we have q^i^ 



3— a „a 



0, A; < z - 1, g" _i = 1, g" = -4 + 2^-°, g; 



i,i+l 



6 



)5— a 



+ 



33-°, and ql^ = {k - i + 2)3^" _ 4(A; - z + l)^-" + 6{k - if-"" - A{k - i - l)^-" + {k 
i — 2)3-°, k > i + 2. For gf^,, k > i + 2, denote m = k — i > 2, then 



9: 



\3-a 



3-a 



\3-a 



i,«+m 



(m + 2)3-° - 4(m + 1)3-° + 6m3-° - 4(m - 1)3-° + (m - 2) 



\3-a 



772 



3-Q 



1 + 1 

m 



3-a 



41 1 + - 
m 



3-Q 



4(1-- 



3-a 



1-1 



3-a 



rri 



3-a 



6 + E 



fc=0 



3 — a 



m 



3-a 



, fc=4 



3 — a 






iy-4fir-4f^V+fz! 

m/ \m/ \ m J \ m 



m' 



a-l 



(3-a)(2-a)(l-a)(-a) 



?Ti^ 



>o, 



noting that the terms (with k being odd numbers) of the series are zero and the terms 
(with k being even numbers) are positive. According to Remark 2.1, all the corresponding 
properties of p°;. are obtained. D 
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